# -*- coding: utf-8 -*-
"""
Created on Mon Sep  6 16:57:38 2021

@author: administer
"""
import numpy as np
import scipy.io as scio
from scipy.fftpack import fft
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

path = 'H:\\研究生\\课题组\\EAPC-WSS'                          #修改路径为你的路径
outpath = 'H:\\研究生\\课题组\\EAPC-WSS\\1'



name = f"/SNR=0-n6-25-0.1-random-pca_0.001.mat"


data = scio.loadmat(path + name)
label = data['label']
y1 = data['data1']
phi1 = data['phi1']

y2 = data['data2']
phi2 = data['phi2']

label_array = np.array(label)

y_array1 = np.array(y1)
phi_array1 = np.array(phi1)

y_array2 = np.array(y2)
phi_array2 = np.array(phi2)

x_data1 = []
x_data2 = []
x_data = []
for i in range(len(label_array)):
    x_data1.append(np.dot(np.linalg.pinv(phi_array1[i]),y_array1[i]))    #np.linalg.pinv意思为计算伪逆
    x_data2.append(np.dot(np.linalg.pinv(phi_array2[i]),y_array2[i]))

x_data1 = np.array(x_data1)
x_data2 = np.array(x_data2)

x_data1_real = np.real(x_data1)
x_data1_imag = np.imag(x_data1)
x_data2_real = np.real(x_data2)
x_data2_imag = np.imag(x_data2)


x_data1_real = x_data1_real.reshape(-1,1)
x_data2_real = x_data2_real.reshape(-1,1)

x_data1_imag = x_data1_imag.reshape(-1,1)
x_data2_imag = x_data2_imag.reshape(-1,1)

x_data3_real = np.concatenate((x_data1_real,x_data2_real),axis = 1)
x_data3_imag = np.concatenate((x_data1_imag,x_data2_imag),axis = 1)



# 选择一个切片进行可视化（假设是第一个样本的频谱数据）
x_real = x_data3_real # 实部
x_imag = x_data3_imag  # 虚部

# 定义时间轴（根据数据情况）
L = 195
R = 1
K = 91
K0 = 10
fnyq = 10e10
TimeResolution = 1 / fnyq
TimeWin = [0, L * R * K - 1, L * R * (K + K0) - 1]
t_axis = np.arange(TimeWin[0], TimeWin[-1] + 1) * TimeResolution

# 定义数字时间轴
Digital_time_axis = np.linspace(t_axis[0], t_axis[-1], x_real.shape[0])

# 要可视化的样本数据（选择一个特定的频率分量，例如第一个）
DigitalSamples1 = x_real
DigitalSignalSamples = x_imag

# 绘制 DigitalSamples1
time_axis_min = 0
time_axis_max = 2e-7

# 振幅范围可以手动设置为较小的范围
amplitude_min = -30000
amplitude_max = 30000

# time_axis_max = Digital_time_axis.max()  # 使用数据的最大值

# 绘图
plt.figure()
plt.plot(Digital_time_axis, DigitalSamples1, 'r')
plt.title('DigitalSamples1')
plt.xlabel('times (s)')
plt.ylabel('magnitude')
plt.xlim(time_axis_min, time_axis_max)
plt.ylim(amplitude_min, amplitude_max)
plt.grid(True)


# 显示图像
plt.show()

# print(x_data3_real.shape)  #118170000 2
pca = PCA(1)
x_data_real = pca.fit_transform(x_data3_real)
x_data_imag = pca.fit_transform(x_data3_imag)

x_data_real = x_data_real.reshape(1*195*101,1)                 #根据你产生的数据集修改大小
x_data_imag = x_data_imag.reshape(1*195*101,1)

for j in range(len(x_data_imag)):
    x_data.append(complex(x_data_real[j,0], x_data_imag[j,0]))
x_data = np.array(x_data)
x_data = x_data.reshape(1,195,101)



x_data_array = np.array(x_data)
x_data_array_fft = fft(x_data_array)


x_data_real = np.real(x_data_array_fft)
x_data_img = np.imag(x_data_array_fft)


x_final = np.stack((x_data_real,x_data_img), axis = 3)

# 选择一个切片进行可视化（假设是第一个样本的频谱数据）
x_real = x_final[0, :, :, 0]  # 实部
x_imag = x_final[0, :, :, 1]  # 虚部

# 定义时间轴（根据数据情况）
L = 195
R = 1
K = 91
K0 = 10
fnyq = 10e10
TimeResolution = 1 / fnyq
TimeWin = [0, L * R * K - 1, L * R * (K + K0) - 1]
t_axis = np.arange(TimeWin[0], TimeWin[-1] + 1) * TimeResolution

# 定义数字时间轴
Digital_time_axis = np.linspace(t_axis[0], t_axis[-1], x_real.shape[0])

# 要可视化的样本数据（选择一个特定的频率分量，例如第一个）
DigitalSamples1 = x_real[:, 0]
DigitalSignalSamples = x_imag[:, 0]

# 绘制 DigitalSamples1
time_axis_min = 0
time_axis_max = 2e-7

# 振幅范围可以手动设置为较小的范围
amplitude_min = -30000
amplitude_max = 30000

# time_axis_max = Digital_time_axis.max()  # 使用数据的最大值

# 绘图
plt.figure()
plt.plot(Digital_time_axis, DigitalSamples1, 'r')
plt.title('DigitalSamples1')
plt.xlabel('times (s)')
plt.ylabel('magnitude')
plt.xlim(time_axis_min, time_axis_max)
plt.ylim(amplitude_min, amplitude_max)
plt.grid(True)


# 显示图像
plt.show()


data_out = {'x':x_final,'label':label}

scio.savemat(outpath + filename,data_out)



